Universaltool Regressionsanalyse II (Mehrebenenmodelle)

Tag 1: Multiple Regression

Samuel Merk

01.04.2022

Meine Gedanken zur Workshopgestaltung

  • Möglichst viel Aktivität bei den Teilnehmer*innen
    • Weniger Vorturnen durch mich, mehr selbst ausprobieren (und erstmal scheitern 😃)
    • Ich zeige zunächst alles in , jeder benutzt danach seine Lieblingssoftware - ich helfe bei Fragen zu JASP/jamovi/SPSS/STATA/etc. so gut ich kann …
  • Möglichst viel Interaktion
    • Fragen zu Begriffen gerne direkt stellen
    • Wünsche wiederholter/weiterer Erklärungen oder Elaborationen gerne direkt äußern
    • Gerne Peer-to-Peer Interaktion in Break-Out-Rooms
    • Gemeinsames Interpretieren der Lehrbuchtexte (in Sitzung 2)
  • Möglichst differenziertes Arbeiten
    • Gestufte Lösungshilfen
    • Break-Out-Rooms für “aktiv Modellierende” und “passiv Nachvollziehende”
    • Kognitiv aktivierende Inputs die auf verschiedenen Verständnisebenen rezipiert werden können
  • Materials

Worked Out Example: Interpretation multipler Regressionsmodelle

Datensatz 1: Kid IQ

Der Datensatz Kid IQ stammt aus einem der empfohlenen Lehrbücher(Gelman & Hill, 2007). Ursprünglich stammt er aus der National Longitudinal Survey of Youth. Wir werden die folgenden Variablen Nutzen:

  • kid_score = Rohwert des Kidnes in einem Intelligenztest
  • mom_hs = Dummyvariable Highschoolabhschluss (1 = Highschoolabschluss, 0 = kein Highschoolabschluss)
  • mom_iq = IQ der Mutter
  • mom_age = Alter der Mutter
  • mom_work =
    • 1: mother did not work in first three years of child’s life
    • 2: mother worked in second or third year of child’s life
    • 3: mother worked part-time in first year of child’s life
    • 4: mother worked full-time in first year of child’s life

Import der Daten

library(tidyverse)
data_kidiq <- read_delim("https://raw.githubusercontent.com/sammerk/did_data/master/kidiq.csv", delim = ";")

Diejenigen, die nicht nutzen wollen, können den Datensatz hier im SPSS-Format herunterladen.

Modelle mit einem metrischen Prädiktor

Parametrisierung

Zunächst soll der kid_score mit der metrischen Variable mom_iq prädiziert werden. In erfolgt das mit der Syntax

mod01 <- lm(kid_score ~ mom_iq, data = data_kidiq)
mod01
## 
## Call:
## lm(formula = kid_score ~ mom_iq, data = data_kidiq)
## 
## Coefficients:
## (Intercept)       mom_iq  
##       25.80         0.61

Grafisch kann dieses Modell wie folgt repräsentiert werden:

library(hrbrthemes)
data_kidiq %>% 
  ggplot(., aes(mom_iq, kid_score)) + 
  geom_point() + 
  stat_smooth(method = "lm", se = F) + 
  theme_ipsum()
## `geom_smooth()` using formula 'y ~ x'

Und die formale Beschreibung lautet:

\[ \operatorname{\widehat{kid\_score}} = 25.8 + 0.61(\operatorname{mom\_iq}) \]

Da der kid_score offensichtlich nicht einer Standardskalierung (z-Werte, IQ-Werte, t-Werte, …) entspricht, kann die Stärke des Effekts nur anhand einer Standardisierung bewertet werden.

mod02 <- lm(scale(kid_score) ~ scale(mom_iq), data = data_kidiq)
mod02
## 
## Call:
## lm(formula = scale(kid_score) ~ scale(mom_iq), data = data_kidiq)
## 
## Coefficients:
##   (Intercept)  scale(mom_iq)  
##    -5.068e-16      4.483e-01

Dies leisten auch packages wie {sjPlot} oder {stargazer} die darüber hinaus auch noch mehrere Modelle vergleichend darstellen können.

library(sjPlot)
tab_model(mod01, show.std = T, show.ci = F, show.p = F)
  kid_score
Predictors Estimates std. Beta
(Intercept) 25.80 -0.00
mom iq 0.61 0.45
Observations 434
R2 / R2 adjusted 0.201 / 0.199

Nach Cohen (1988) gilt dabei:

  • \(\beta \approx .1 \Rightarrow \text{kleiner Effekt}\)
  • \(\beta \approx .3 \Rightarrow \text{moderater Effekt}\)
  • \(\beta \approx .5 \Rightarrow \text{starker Effekt}\)

Inferenzstatistik

Unser \(\beta = .45\) beschreibt lediglich die Steigung einer Geraden, wenn man diese nach dem Kriterium der kleinsten Quadrate auf dem Datensatz kid_iq optimiert. Möchte man Schlussfolgerungen über den datengenerierenden Mechanismus anstellen benötigt man Inferenzstatistik.
p-Werte kann man etwa anhand der summary()-Funktion oder der tab_model() Funktion erhalten.

summary(mod01)
## 
## Call:
## lm(formula = kid_score ~ mom_iq, data = data_kidiq)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -56.753 -12.074   2.217  11.710  47.691 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 25.79978    5.91741    4.36 1.63e-05 ***
## mom_iq       0.60997    0.05852   10.42  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 18.27 on 432 degrees of freedom
## Multiple R-squared:  0.201,  Adjusted R-squared:  0.1991 
## F-statistic: 108.6 on 1 and 432 DF,  p-value: < 2.2e-16

Dabei sind zwei Inferenzstatistiken enthalten:

  1. Jeder Koeffizient \(\beta_i\) wird gegen die 0 getestet.
  2. Das Gesamtmodell wird gegen \(R^2 = 0\) getestet.

JZS-Bayes-Faktoren sind via {BayesFactor} ermittelbar:

library(BayesFactor)
lmBF(formula = kid_score ~ mom_iq, data = data_kidiq)
## Bayes factor analysis
## --------------
## [1] mom_iq : 4.516119e+19 ±0%
## 
## Against denominator:
##   Intercept only 
## ---
## Bayes factor type: BFlinearModel, JZS

Diagnostik der Voraussetzungen für die Inferenzstatistik

Diese Inferenzstatistiken sind nur unter Annahmen über den datengenerierenden Mechanismus berechenbar. Sind diese Annahmen verletzt sind die Inferenzstatistiken (mehr oder weniger stark) verzerrt. Daher nimmt die Diagnostik der Annahmen eine zentrale Bedeutung ein. Angenommen werden muss (mindestens)

  • Linearität der Beziehung
  • Normalverteilung der Residuen
  • Homoskedastizität der Residuen
  • Unabhängigkeit der Residuen

Eine gute Grundlage für die Diagnostik bietet das {performance}-package:

library(performance)
check_model(mod01)

Interpretation

Regressionsmodelle beschreiben immer nur bedingte Erwartungswerte von Datenpunkten (“für eine Gruppe von Müttern deren IQ im Schnitt X ist, ist ein kid_score von X am wahrscheinlichsten”). Woher diese Wahrscheinlichkeit rührt und wie gut mit dieser eine kausale Relationierung von X und Y gerechtfertigt werden kann hängt maßgeblich vom Design der Studie ab!

Modelle mit einem dichotomen Prädiktor (Dummyvariable)

Parametrisierung

Die Regression ist insofern ein recht universelles Modellierungstool, als dass es auch die Aufnahme dichotomer Prädiktoren erlaubt. Bspw. kann man mit der Variable mom_hs der kid_score prädizieren:

ggplot(data_kidiq, aes(mom_hs, kid_score)) + 
  geom_point(alpha = .3) + 
  stat_smooth(method = "lm", se = F) + 
  theme_ipsum()
## `geom_smooth()` using formula 'y ~ x'

Die resultierenden Koeffizienten stellen die arithmetischen Mittelwerte der Gruppen dar und der p-Wert von \(\beta_1\) tested die \(H_0: EW(Gruppe_1) = EW(Gruppe_2)\).

mod03 <- lm(kid_score ~ mom_hs, data = data_kidiq)
mod03
## 
## Call:
## lm(formula = kid_score ~ mom_hs, data = data_kidiq)
## 
## Coefficients:
## (Intercept)       mom_hs  
##       77.55        11.77

Wobei \[ \operatorname{kid\_score} = \beta_{0} + \beta_{1}(\operatorname{mom\_hs}) + \epsilon \]

Diagnostik der Voraussetzungen für die Inferenzstatistik

Da Regressionsmodelle mit Dummyvariablen eben auch Regressionsmodelle sind 😃 gilt es auch bei diesen dieselben Voraussetzungen zu prüfen:

check_model(mod03)

Modelle mit mehreren Prädiktoren (multiple Regression)

Parametrisierung

In Regressionsmodelle können problemlos mehrere Prädiktoren aufgenommen werden.

mod04 <- lm(kid_score ~ mom_iq + mom_hs, data = data_kidiq)
summary(mod04)
## 
## Call:
## lm(formula = kid_score ~ mom_iq + mom_hs, data = data_kidiq)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -52.873 -12.663   2.404  11.356  49.545 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 25.73154    5.87521   4.380 1.49e-05 ***
## mom_iq       0.56391    0.06057   9.309  < 2e-16 ***
## mom_hs       5.95012    2.21181   2.690  0.00742 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 18.14 on 431 degrees of freedom
## Multiple R-squared:  0.2141, Adjusted R-squared:  0.2105 
## F-statistic: 58.72 on 2 and 431 DF,  p-value: < 2.2e-16

Grafisch kann dies dann im entsprechenden n-dimensionalen Raum repräsentiert werden (siehe Folien aus Universaltool Regressionsanalyse I). In der formalen Repräsentation wird schlicht das Polynom um einen Summanden erweitert: \[ \operatorname{kid\_score} = \beta_{0} + \beta_{1}(\operatorname{mom\_iq}) + \beta_{2}(\operatorname{mom\_hs}) + \epsilon \]

Interpretation

Besonders interessant an multiplen Regressionsmodellen ist, dass bestimmte kausale Mechanismen falsifiziert werden können. Dazu werden dann meist die Ergebnisse der einfachen Modelle (mit jeweils einem Prädiktor) mit den Parameterschätzungen im multiplen Modell verglichen:

tab_model(mod01, mod03, mod04, show.ci = F)
  kid_score kid_score kid_score
Predictors Estimates p Estimates p Estimates p
(Intercept) 25.80 <0.001 77.55 <0.001 25.73 <0.001
mom iq 0.61 <0.001 0.56 <0.001
mom hs 11.77 <0.001 5.95 0.007
Observations 434 434 434
R2 / R2 adjusted 0.201 / 0.199 0.056 / 0.054 0.214 / 0.210

So können Annahmen über Supressor- oder Konfounderkonstellationen falsifiziert werden. Hier etwa wird klar, dass die Daten mit dem Mechanismus A nicht im Einklang stehen, aber mit dem Mechanismus B und C.

Diagnostik der Voraussetzungen für die Inferenzstatistik

Auch Modelle der multiplen Regression unterliegen den gleichen 4 Annahmen. Hinzu kommt (je nach Methode der Schätzung), dass keine Kolinearität vorliegen sollte. check_model() nimmt dies Annahme direkt in die Prüfung mit auf, wenn das Argument der Funktion ein Modell mit mehreren Prädiktoren darstellt.

check_model(mod04)

Typischerweise werden dann wieder p-Werte für die Hypothesen \(\beta_i = 0\) und $ R^2 = 0$ berechnet.
Bayes Faktoren beziehen sich jedoch immer auf das inkrementelle \(R^2\), also die Frage: Klärt ein Prädiktor oder eine Gruppe von Prädiktoren zuätzliche Varianz auf?

# Test gegen intercept only modell
lmBF(formula = kid_score ~ mom_iq + mom_hs, data = data_kidiq)
## Warning: data coerced from tibble to data frame
## Bayes factor analysis
## --------------
## [1] mom_iq + mom_hs : 1.508721e+20 ±0%
## 
## Against denominator:
##   Intercept only 
## ---
## Bayes factor type: BFlinearModel, JZS
# Test gegen Modell mit mom_hs als Prädiktor
lmBF(formula = kid_score ~ mom_iq + mom_hs, data = data_kidiq)/
  lmBF(formula = kid_score ~ mom_hs, data = data_kidiq)
## Warning: data coerced from tibble to data frame

## Warning: data coerced from tibble to data frame
## Bayes factor analysis
## --------------
## [1] mom_iq + mom_hs : 7.578673e+15 ±0%
## 
## Against denominator:
##   kid_score ~ mom_hs 
## ---
## Bayes factor type: BFlinearModel, JZS
# Test gegen Modell mit mom_iq als Prädiktor
lmBF(formula = kid_score ~ mom_iq + mom_hs, data = data_kidiq)/
  lmBF(formula = kid_score ~ mom_iq, data = data_kidiq)
## Warning: data coerced from tibble to data frame

## Warning: data coerced from tibble to data frame
## Bayes factor analysis
## --------------
## [1] mom_iq + mom_hs : 3.340747 ±0%
## 
## Against denominator:
##   kid_score ~ mom_iq 
## ---
## Bayes factor type: BFlinearModel, JZS

Modelle mit Interaktionseffekt eines metrischen und eines dichotomen Prädiktors

Parametrisierung

Interessiert man sich dafür, wie sich die Effekte je nach Ausprägung einer weiteren nominalen Variablen unterscheiden, kann man Modelle mit Interaktionen eines metrischen und eines dichotomen Prädiktors spezifizieren.

ggplot(data_kidiq, aes(mom_iq, kid_score, 
                       color = as.factor(mom_hs), 
                       group = as.factor(mom_hs))) + 
  geom_point() + 
  stat_smooth(method = "lm") + 
  theme_ipsum()
## `geom_smooth()` using formula 'y ~ x'

Die entsprechenden Koeffizienten des Modells mit Interaktionseffekt sind jedoch wesentlich schwieriger zu interpretieren:

mod05 <- lm(kid_score ~ mom_hs + mom_iq + mom_hs:mom_iq, data = data_kidiq)
mod05
## 
## Call:
## lm(formula = kid_score ~ mom_hs + mom_iq + mom_hs:mom_iq, data = data_kidiq)
## 
## Coefficients:
##   (Intercept)         mom_hs         mom_iq  mom_hs:mom_iq  
##      -11.4820        51.2682         0.9689        -0.4843

\[ \operatorname{kid\_score} = \alpha + \beta_{1}(\operatorname{mom\_hs}) + \beta_{2}(\operatorname{mom\_iq}) + \beta_{3}(\operatorname{mom\_hs} \times \operatorname{mom\_iq}) + \epsilon \]

Leichter wird dies, wenn man den kontinuierlichen Prädiktor zentriert oder standardisiert.

library(equatiomatic)
data_kidiq <- data_kidiq %>% 
  mutate(mom_iq_centered = as.numeric(scale(mom_iq, center = T, scale = F)),
         mom_iq_zstand = as.numeric(scale(mom_iq, center = T, scale = T)),
         mom_age_zstand = as.numeric(scale(mom_age, center = T, scale = T)))

mod06 <- lm(kid_score ~ mom_hs + mom_iq_zstand + mom_hs:mom_iq_zstand, data = data_kidiq)

mod06
## 
## Call:
## lm(formula = kid_score ~ mom_hs + mom_iq_zstand + mom_hs:mom_iq_zstand, 
##     data = data_kidiq)
## 
## Coefficients:
##          (Intercept)                mom_hs         mom_iq_zstand  
##               85.407                 2.841                14.533  
## mom_hs:mom_iq_zstand  
##               -7.264
ggplot(data_kidiq, aes(mom_iq_zstand, kid_score, 
                       color = as.factor(mom_hs), 
                       group = as.factor(mom_hs))) + 
  geom_point() + 
  stat_smooth(method = "lm") + 
  theme_ipsum()
## `geom_smooth()` using formula 'y ~ x'

extract_eq(mod06, use_coefs = T, terms_per_line = 2, wrap = T)

\[ \begin{aligned} \operatorname{\widehat{kid\_score}} &= 85.41 + 2.84(\operatorname{mom\_hs})\ + \\ &\quad 14.53(\operatorname{mom\_iq\_zstand}) - 7.26(\operatorname{mom\_hs} \times \operatorname{mom\_iq\_zstand}) \end{aligned} \]

Diagnostik der Voraussetzungen für die Inferenzstatistik

Auch hier gelten dieselben Annahmen, die wieder gut mit check_model() evaluierbar sind:

check_model(mod06)

Modelle mit Interaktionseffekt zweier metrischer Prädiktoren

Betrachtet man Modelle mit Interaktionen zwischen zwischen zwei metrischen Variablen, sind diese am besten interpretierbar, wenn beide zentriert sind.

mod07 <- lm(kid_score ~ mom_iq_zstand*mom_age_zstand, data = data_kidiq)
mod07
## 
## Call:
## lm(formula = kid_score ~ mom_iq_zstand * mom_age_zstand, data = data_kidiq)
## 
## Coefficients:
##                  (Intercept)                 mom_iq_zstand  
##                      86.8761                        9.2097  
##               mom_age_zstand  mom_iq_zstand:mom_age_zstand  
##                       1.1404                       -0.8633
extract_eq(mod07, use_coefs = T, wrap = T)

\[ \begin{aligned} \operatorname{\widehat{kid\_score}} &= 86.88 + 9.21(\operatorname{mom\_iq\_zstand}) + 1.14(\operatorname{mom\_age\_zstand}) - 0.86(\operatorname{mom\_iq\_zstand} \times \operatorname{mom\_age\_zstand}) \end{aligned} \] In diesen Modellen können dann die sogenannten Haupteffekte, also die \(\beta_i\) der einfachen Summanden ohne multiplikative Terme (hier \(\beta_1\) und \(\beta_2\)) als die Effekte der Variable \(i\) gesehen werden, wenn die andere Variable eine durchschnittliche Ausprägung (nach Zentrierung = 0) hat.

Richtig hübsch geplottet wird das per default mit der Funktion plot_model()

plot_model(mod07, type = "int") +
  theme_ipsum()

Übung: Multiple Regression mit Dummyvariablen

Datensatz 2: Effekte der Klassengröße und Wertzuschreibung

Das Student Teacher Achievement Ratio (STAR) Projekt ist eines der größten Experimente zu Effekten der Klassengrößenreduktion auf die Schüler*innenleistung. Die Klassen wurden randomisiert drei Bedingungen zugewiesen: Große Klasse, kleine Klasse, große Klasse mit Hilfslehrkraft. Für die dritte Klasse findet sich diese Information in der Variable g3classtype (1 = SMALL CLASS, 2 = REGULAR CLASS, 3 = REGULAR + AIDE CLASS). Außerdem wurde über einen Lehrerfragebogen die Wertzuschreibung der Schülerinnen und Schüler erfasst. Im Datensatz ist die Variable mit dem Kürzel g4ptvalu zu finden. Ein Beispielitem lautet This student thinks that school is important mit den Antwortmöglichkeiten 1 = Never, 2, 3 = Sometimes, 4, 5 = Always. Die Mathematikleistung ist für die dritte und vierte Klasse in den Variablen g3tmathss und g4tmathss enkodiert.

Die Daten können wie folgt heruntergeladen werden (ein .sav-fiel gibt es hier)

data_star <- read_spss("https://raw.githubusercontent.com/sammerk/did_data/master/STAR.sav") 

Fragestellung 1

Mein Vorschlage wäre zunächst die folgende Fragestellung zu bearbeiten (Lösungshinweise und die Musterlösung sind dann jeweils hinter den Tabsets zu finden)

Fragestellung

Welchen Effekt zeigt die Klassengrößenreduktion auf die Mathematikleistung? Schätzen Sie die Effekte, evaluieren Sie Effektstärken und Inferenzstatistiken und überlegen inwiefern die Effekte kausal interpretiert werden können.

Modellierung

Hinweis 1

Doe Klassengröße g4classsize stellt ja einen nominalen Prädiktor mit 3 Ausprägungen dar. Daher macht ein Regressionsmodell mit zwei Dummyvariablen sinnd

Lösung

mod08 <- lm(g3tmathss ~ as.factor(g3classtype), data = data_star)
summary(mod08)
## 
## Call:
## lm(formula = g3tmathss ~ as.factor(g3classtype), data = data_star)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -135.958  -27.275   -1.275   25.921  158.725 
## 
## Coefficients:
##                         Estimate Std. Error t value Pr(>|t|)    
## (Intercept)             622.9577     0.9017 690.845  < 2e-16 ***
## as.factor(g3classtype)2  -6.8787     1.2928  -5.321 1.07e-07 ***
## as.factor(g3classtype)3  -7.6824     1.2237  -6.278 3.66e-10 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 39.7 on 6074 degrees of freedom
##   (5524 Beobachtungen als fehlend gelöscht)
## Multiple R-squared:  0.007415,   Adjusted R-squared:  0.007088 
## F-statistic: 22.69 on 2 and 6074 DF,  p-value: 1.528e-10

Modellinterpretation

Das Intercept stellt das arithmetische Mittel in der Referenzgruppe - also den kleinen Klassen - dar. Die Slopes die Abweichungen der Mittelwerte der anderen Gruppen von diesen Mittelwerten. Große Klassen und große klassen mit Hilfslehrkraft performen also signifikant schlechter, wobei die Effekte nach Cohen klein sind und dekriptiv kein Unterschied zwischen den beiden großen Klassentypen zu sehen ist.

Modelldiagnostik

check_model(mod08)

Normality sieht ganz gut aus. Die großen Abstände und CI bei Linearity und Homogenity rphren von der Dummykodierung, sind aber auch unproblematisch.
Problematisch ist aber der Mehrebenenkontext (Schüler*innen in Klassen etc.). Dies führt zu einer Abhängigkeit der Residuen welche man etwa mit dem ICC1 quantifizieren kann (siehe nächste Sitzung).

psychometric::ICC1.lme(g3tmathss, g3schid, data = data_star)
## [1] 0.1847848

Fragestellung 2

Fragestellung

Welchen prädiktiven Effekt zeigt die Wertzuschreibung in Klasse 4 g4ptvalu auf die Mathematikleistung in Klasse 4 (Variable g4tmathss unter Kontrolle der Mathematikleistung in Klasse 3 (Variable g3tmathss)? Schätzen Sie die Effekte, evaluieren Sie Effektstärken und Inferenzstatistiken und überlegen Sie welche hypothetischen Kausalmechanismen Sie mit den Daten widerlegen/belegen können.

Modellierung

Hinweis 1

Die Formulierung “unter Kontrolle von” ist typisch für multiple Regressionsmodelle. Um potentiell Confounding oder Suppression Konstellationen finden zu können, macht es Sinn, erst einfache Regressionen mit den Prädiktoren zu schätzen und dann ein Modell, das beide Prädiktoren enthält.

Lösung

mod09 <- lm(g4tmathss ~ g3tmathss, data = data_star)
mod10 <- lm(g4tmathss ~ g4ptvalu, data = data_star)
mod11 <- lm(g4tmathss ~ g3tmathss + g4ptvalu, data = data_star)

tab_model(mod09, mod10, mod11, show.ci = F, show.std = T)
  TOTAL MATH SCALE SCORE
CTBS GRADE 4
TOTAL MATH SCALE SCORE
CTBS GRADE 4
TOTAL MATH SCALE SCORE
CTBS GRADE 4
Predictors Estimates std. Beta p Estimates std. Beta p Estimates std. Beta p
(Intercept) 209.29 0.00 <0.001 613.36 -0.00 <0.001 230.39 -0.00 <0.001
TOTAL MATH SCALE SCORE
SAT GRADE 3
0.80 0.70 <0.001 0.70 0.64 <0.001
GRADE 4 PARTICIPATION
SUBSCORE: VALUE
7.49 0.36 <0.001 3.20 0.16 <0.001
Observations 4069 1850 1793
R2 / R2 adjusted 0.483 / 0.483 0.132 / 0.131 0.490 / 0.489

Modellinterpretation

Der Effekt der Wertzuschreibung sinkt beträchtlich unter Hinzunahme der Leistung in Klasse drei, wohingegen die prädiktive Kraft der Leistung in Klasse 3 unter Hinzunahme der Wertvariable kaum sinkt.
Dies ist im Einklang mit der Annahme, dass der Effekt der Wertvariable durch die Leistungs in Klasse vier konfundiert ist.

set.seed(12345)
D <- dagify(L4 ~ W4 + L3,
            W4 ~ L3)

ggdag(D) + 
  theme_dag_blank() + 
  ggtitle("Konfundierung des Effekt der Wertzuschreibung auf die Leistung in Klasse 4", "durch die Lesitung in Klasse 3")

Jedoch stehen die Daten auch mit einem Modell in Einklang, indem die Wertzuschreibung und die Leistung von einer weiteren unbeobachteten Variable beeinflusst sind.

set.seed(12345)
E <- dagify(L4 ~ W4 + L3,
            W4 ~ U,
            L3 ~ U)

ggdag(E) + 
  theme_dag_blank() + 
  ggtitle("Wertzuschreibung und die Leistung in Klasse 4", "sind von einer weiteren unbeobachteten Variable beeinflusst")

Modelldiagnostik

check_model(mod11)

Auch hier ist die Mehrebenenkonstellation wieder das zentrale Problem und leider nicht in den Plots sichtbar.

Übung: Multiple Regression mit Interaktionseffekten

Datensatz 3: European study on family care of older people (efc)

Aus dem Datensatz der European study on family care of older people wollen wir die Variablen c12hour barthtot und c161sex betrachten. neg_c_7 misst den negativen Impact der Pflegeaufgabe mit einer Skala, die aus 7 Likert Items (z.B. Does caregiving have a negative effect on your emotional well-being?, Balducci et al., 2008) besteht. c12hour stellt die durchschnittliche Anzahl Pflegearbeit pro Woche dar und barthtot den sogenannten Barthel Index. Dieser beschreibt inwiefern Aktivitäten des täglichen Lebens eingeschränkt sind (größere Werte = größere Einschränkungen). c161sex stellt eine dichotome Gendererfassung (1 = Male, 2 = Female) dar.

Die Daten können wie folgt heruntergeladen werden (ein .sav-file gibt es hier)

data_efc <- read_spss("https://raw.githubusercontent.com/sammerk/did_data/master/data_efc.sav") 

Fragestellung 1

Fragestellung

Gibt es differentielle Effekte des Pflegeumfangs auf den negativen Impact je Geschlecht?

Modellierung

Hinweis 1

Die Formulierung “differentieller Effekt” kann grafisch als unterschiedlich stark steigende Regressionsgeraden je Geschlecht interpretiert werden. Dies erzielt man durch multiplikation der beiden Prädiktoren.

Hinweis 2

Wichtig kann dabei sein, die Variablen zu rekodieren. c161sex ist bspw. mit 1 und 2 kodiert und c12hour hat eine natürliche Metrik.

Lösung

mod12 <- lm(neg_c_7 ~ scale(c12hour)*as.factor(c161sex), data = efc)
summary(mod12)
## 
## Call:
## lm(formula = neg_c_7 ~ scale(c12hour) * as.factor(c161sex), data = efc)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -7.3724 -2.7187 -0.7549  1.7168 16.2755 
## 
## Coefficients:
##                                    Estimate Std. Error t value Pr(>|t|)    
## (Intercept)                        11.52122    0.25753  44.737  < 2e-16 ***
## scale(c12hour)                      1.06961    0.26342   4.060 5.33e-05 ***
## as.factor(c161sex)2                 0.44121    0.29467   1.497    0.135    
## scale(c12hour):as.factor(c161sex)2 -0.09479    0.29966  -0.316    0.752    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 3.73 on 887 degrees of freedom
##   (17 Beobachtungen als fehlend gelöscht)
## Multiple R-squared:  0.07012,    Adjusted R-squared:  0.06697 
## F-statistic: 22.29 on 3 and 887 DF,  p-value: 6.36e-14

Modellinterpretation

Während c12hour einen signifikanten Effekt moderater Größe zeigt, ist dies für die Gendervariable und deren Interaktion nicht der Fall. Aber Achtung: Dies ist keine Evidenz für die Abwesenheit dieser Effekte. Selbige kann man durch folgende BF erhalten:

lmBF(formula = neg_c_7 ~ c12hour_zstand*c161sex_I,
     data = efc %>% 
       select(neg_c_7, c12hour, c161sex) %>% 
       mutate(c161sex_I = c161sex - 1,
              c12hour_zstand = as.numeric(scale(c12hour))) %>% 
       na.omit()) /
  lmBF(formula = neg_c_7 ~ c12hour_zstand,
     data = efc %>% 
       select(neg_c_7, c12hour, c161sex) %>% 
       mutate(c161sex_I = c161sex - 1,
              c12hour_zstand = as.numeric(scale(c12hour))) %>% 
       na.omit()) 
## Bayes factor analysis
## --------------
## [1] c12hour_zstand * c161sex_I : 0.03676294 ±0%
## 
## Against denominator:
##   neg_c_7 ~ c12hour_zstand 
## ---
## Bayes factor type: BFlinearModel, JZS

Die Funktion lmBF() kann allerdings nicht scale() und as.factor() als Argumente enthalten. Daher müssen diese zuvor rekodiert werden. Ebenfalls hat lmBF() keine “silent case wise deletion”, weshalb die Missings händisch gelöscht werden müssen.

Modelldiagnostik

check_model(mod12)

Fragestellung 2

Fragestellung

Moderiert der Barthel Index den Effekt des Pflegeumfangs auf den negativen Impact?

Modellierung

Hinweis

“Moderiert”, “interagiert”, “differiert ein Effekt” etc. sind typische Ausdrucksweisen für die Interpretation multiplikativer Terme in Regressionsgleichungen.

Lösung

mod13 <- lm(neg_c_7 ~ c12hour * barthtot, data = efc)
tab_model(mod13, show.ci = F, show.std = T)
  Negative impact with 7
items
Predictors Estimates std. Beta p std. p
(Intercept) 15.56 0.05 <0.001 0.177
average number of hours
of care per week
-0.01 0.12 0.168 0.002
Total score BARTHEL INDEX -0.06 -0.39 <0.001 <0.001
c12hour:barthtot 0.00 0.09 0.002 0.002
Observations 874
R2 / R2 adjusted 0.180 / 0.178

Modellinterpretation

Alle drei Effekte sind signifikant. Die Größe insbesondere des Interaktionseffekt ist wesentlich leichter zu interpretieren wenn man die Daten plottet oder standardisiert:

plot_model(mod13, type = "pred",
           terms = c("c12hour", 
                     # `barthtot [30,50,70]` definiert **welche** differentiellen 
                     # Effekte dargestellt werden
                     "barthtot [30,50,70]")) + 
  theme_ipsum()

Modelldiagnostik

check_model(mod13)

Übung: Einfache Regression mit Count Data

Datensatz 4: Demand for medical care

Deb & Trivedi (1997) befragten über 4000 Individuen die älter als 60 waren (u.a.) zur number of chronic conditions numchron und der number of physicians office visits ofp, welche als Proxi des demand for medical care diente.
Die Daten können wie folgt heruntergeladen werden (ein .sav-file gibt es hier)

data_medcaredemand <- read_spss("https://raw.githubusercontent.com/sammerk/did_data/master/data_medcaredemand.sav") 

cor(data_medcaredemand$numchron, data_medcaredemand$ofp, method = "spearman")

Fragestellung 1

Fragestellung

Wie stark sind die beiden Variablen numchron und ofp assoziiert?

Lösung

mod14 <- lm(ofp ~ numchron, data = data_medcaredemand)
tab_model(mod14, show.std = T)
  ofp
Predictors Estimates std. Beta CI standardized CI p
(Intercept) 3.75 -0.00 3.46 – 4.04 -0.03 – 0.03 <0.001
numchron 1.31 0.26 1.17 – 1.45 0.23 – 0.29 <0.001
Observations 4406
R2 / R2 adjusted 0.069 / 0.068

Modellinterpretation

Wir sehen einen kleinen bis moderaten signifikanten Effekt. Die kausale Relationierung der Variablen bleibt allerdings völlig ungeklärt.

Modelldiagnostik

check_model(mod14)

Wieder führt die diskrete Prädiktorvariable zu Artefakten in der Referenceline in den ersten beiden Abbildungen. Weit problematischer Ist aber die Abweichung der Verteilung der Residuen von der Normalverteilung. Diese ist typisch für Count Data welche auch schon theoretisch die Linearitätsannahme in Frage stellen. Glücklicherweise gibt es dafür spezielle Regressionsmodelle - z.B. die Poissonregression.

Literatur

Balducci, C., Mnich, E., McKee, K. J., Lamura, G., Beckmann, A., Krevers, B., Wojszel, Z. B., Nolan, M., Prouskas, C., Bien, B., & Oberg, B. (2008). Negative Impact and Positive Value in Caregiving: Validation of the COPE Index in a Six-Country Sample of Carers. The Gerontologist, 48(3), 276–286. https://doi.org/10.1093/geront/48.3.276
Cohen, J. (1988). Statistical power analysis for the behavioral sciences (2. Aufl.). Lawrence Erlbaum.
Deb, P., & Trivedi, P. K. (1997). Demand for medical care by the elderly: A finite mixture approach. Journal of applied Econometrics, 12(3), 313–336. Gelman, A., & Hill, J. (2007). Data analysis using regression and multilevel/hierarchical models. Cambridge University Press.
Nye, B., Hedges, L. V., & Konstantopoulos, S. (1999). The long-term effects of small classes: A five-year follow-up of the Tennessee class size experiment. Educational Evaluation and Policy Analysis, 21(2), 127–142. https://doi.org/10.3102/01623737021002127